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Abstract 


A model for the recovery of gravity anomalies from high precision 
altimeter data is derived which consists of small correction terms to the inverse 
Stokes’ formula. The infinence of ?mknown sea surface topography in the case of 
meandering currents such as theGulf stream is discussed, A formula was derived 
in order to estimate the accuracy of the gravity anomalies from the known accuracy 
of the altimeter data. It is shown that for the case of known harmonic coefficients 
of lower order the range of integration in Stokes inverse formula can be reduced 
very much. 
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1. Introduction 


As part of its Earth and Ocean Physics Application Program (EOPAP) 
the National Aeronautics and Space Administration (NASA) plans in the next ten 
years the launch of some satellites equipped with altimeter for ranging to the 
ocean surface. The announced accuracy of the future altimeter systems lies in 
the range of 10 cm, 

Altimeter measurements of this accuracy will give an important increase 
in the knowledge of the features of the gravity potential of the earth, because the 
ocean surface is not far apart from an equipotential surface of the gravity potential. 
In the present report a model for the recovery of gravity anomalies from this high 
precision altimeter data will be developed. 

ha a study of this kind three questions should be answered: 


1) What is the mathematical model for the solution of the 
problem, which is reasonably simple but sufficiently 
accurate in view of tho given accuracy of tho data? 

2) Can this theoretical model be transformed into a 
stable numerical procedure, suited for a modern 
computer? 

3) What is the accuracy of the result that can he expected 
in view of the anntjunced accuracy of the data ? 


The next section is dedicated to a detailed description of the theoretical 
model. It is based on a solution of Stokes problem, which was developed in a jjre- 
paratory study (Lelgemann, 1975). The present procedure uses the formulae, 
which were derived there, in an inverse application. With this idea the solution 
consists of the three steps: 


a) Correct the altimeter measurements by some small 


1 - 



terms to a suited quantity Na , 


( 1 - 1 ) 



y, . . . , . normal gravity at equator 


b) Use the well-known inverse Stokes' formula for a 
spherical problem 


( 1 - 2 ) 





4Tta 


j 


- l^p)dff ) 


a semi major axis 


c) Correct the result Ag* by some small terms to the 
desired gravity anomaly Ag at the earth's surface. 


A detailed explanation of the method and the terms in formula (1-1) 
and (1-2) is given in the next section. The method is well suited to the re- 
duction of our complicated problem to a solution of a spherical problem. 

As there are no numerical problems in the computation of the correction 
terms the answer to question two is mostly concerned with the spherical solution. 

In the case of spherical approximation much work has already been done in answer- 
ing this question. A direct and also indirect solution of formula (1-2) is discussed 
by (Gopalaplllai, 1974)„ (Smith, 1974) and (Rapp, 1974) both applied the method 
of least squares estimation to the spherical problem. 

The main work done in the present study is the consideration of the 
accuracy, which can be expected under some favourable assumptions about the 
size of the measuring error and the distribution of the measurements. Under 
favourable circumstances ( a good Itnowledge of the satellite orbit and so on) there 
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remain in any case ihree error sources; 


a) Influence of the unlmown sea surface topography. We 
shall discuss in section three the size and influence of 
this error type, 

b) The deficiency of suitable altimeter data on the con- 
tinents. We may overcome this difficulty by neglect- 
ing the outer zones in the integral ( 1 - 2 ); that means, 
we will integrate only up to a certain distance 

For this reason we have to diseuss the error due to 
the neglected parts of the potential function itself. 

c) Influence of measuring errors. We wUl assume that 
mean values in b (k m) * b (k m) blocks are given 
with an accuracy of m(N), e. g. Ft, = b * b - 

10 km * 10 km and m(N) - ilOcm, 


The problem under consideration is the behavior of these errors when 
they are transformed by a linear integra!. operator over the sphere, especially 
in view of the feet, that the kernel 


(1-3) 


(2 sin 0/2)^ 

ip spherical distance 

io • . • • ■ distance on the unit sphere 


of the integral is singular, so that this integral exists only under some restric- 
tions. 


Problems of this form are discussed within statistical methods in 
physical geodesy (Heiskanen-Moritz, 1967). A deep examination of the problem 
of the error propagation in the case of integral operators over the unit sphere 
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was given by Meissl (Meissl, 1971). From both these publicatioi^ we will 
recall some basic relations and definitions about functions on the sphere, which 
we need in our special analysis, 

A function on the surface of the sphere, which is at least Lebesgue-in- 
tegrable, can be e:q)anded into a convergent series of spherical Imrmonics by 

(1-4) f(0,X) 

&=0 

Rojj(0, X), Sna(0, X) fully normalized spherical harmonics. 


The terms 

n 

(1-5) 

W=0 




[aiajRji!a(0, X) + ^nnsSnnC^J 3 


are called the degree variances of the function f, 


(1-6) 


cov(f, £) 



the covariance function of the function f (Heiskanen- Moritz, 1967, (7-19) ) and 


(1-7) var(f^) = II f II® = ^ ftt® 

n = 3 


the variance (or the norm square) of the fiinction f. Because it gives all infor- 
mation about the statistical behavior of the disturbing potential, the covariance 
function cov (N, N) or cov (Ag, Ag) is a basic tool of the error considerations in 
physical geodesy. 

At.present the best information about the statistical behavior of the dis- 
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turbing potential on the geoid is probably the covariance function given by 
Tscherning and Rapp (Tscherning-Rapp, 1974), 


(1-8) 


with 


cov (N, N) = ^ 

n=3 



A 

(n - 1) (n - 2) (n + B) 


Pn(COSt!)), 


A = 425. 28 mgal® 
B - 24 
S = 0.999617 . 


A covariance function of this form presupposed that the function N 
is analytic on the sphere (that is on the geoid). in the aeries for the covariance 
function, this presupposition is expressed by the factor (st”^^)), which runs very 
fest to zero. 

We will try to malce some specifications about the covariance function 
of the error function u. It is well known that in the case of uncorrelated point 
measurements the error covariance function has a certain value at the origin 
and is zero elsewhere. In accordance with our presupposition we have as 
measurements N mean values in certain blocks of area F^, that is, the func- 
tion N is giver, as a step hmetion with discontinuities at the boundaries of the 
blocks. In this case we have to consider not individual measuring errors rather 
an error function 


(1-9) 


p N - N. 


This error function 1/ must be lii^ewise a bounded step ftinctiou. Therefore u 
is Lebesgue-integrable over the sphere with the series evaluation 


V (0j X) — ^ 2 [OflBRna{®» ■*' flnn^HB(®» ^)1 


n cuzO 


-ii- 


( 1 - 10 ) 



under the condition 


( 1 - 11 ) 


m®(N) 


11 ^ 



CO 


9 


where is defined as 

(1"12) ^ (0^® + Kifl) . 


We need both the covariance function of the potential function itself and the 
error covariance function in our discussion about the error propagation. 

V^ithln the estimate of the error propagation in the case of an integral 
operator on the sphere, a very useful tool is the conception of the greatest lower 
bound and the least upper bound (Meissl, 1971, p. 26). For example, wc can 
Write for the norm between N and Ag 


(1-13) 


GLB( ]\n|) < 



LUB( |Xnl). 


GLB is the abbreviation of greatest lower bound, some sort of a generalized 
minimum and LUB stands for least upper bound, which is some sort of ge-, 
neralized maximum. Under very general conditions the GLB and LUB of 
the norm ratios of the main functions of gravimetric geodesy are computed 
by Meissl (Meissl, 1971). 

In our Special problem we have 


(1-14) 


N 

-i-'ino 


2l. . 1 

78 (n - 1) 
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Since we know the covariance function of the disturbing potential itself we can 
compute the norm ratio explicitly by 


(1-15) 



Using the covariance fimction given by (Tscheriung-Rapp, 1974) we get 


(1-16) 


a 




because the disturbing potential is approximated by an analytical function. 

In the case of the error function we have, corresponding to (1-14), 
the relationship 


(1-17) 



e error of the gravity anomalies Ag 


with an estimate of the norm ratio of 


(1-18) 



(see (Meissl, 1971, page 46) ). The difference in the expressions (1-16) and 
(1-18) is of course due to the fact, that the error function v is not an analytical 
function. 


We see from (1-18), that small errors In the data may accumulate in 
unfavorable cases to infinite large errors in the wanted result. There are two 
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principal possibilities to overcome this undesirable behaviort 


a) Smoothing out also the error fimction V to an analyti- 
cal function. This is possible by a suitable interpolation 
procedure. 

b) Smoothing out the operator. 


jVccording to the assumed form of our measurements as mean values of certain 
blocks we shall use in our accuracy considerations the second possibility. M 
geodetic problems this is usually done by a separate treatment of a certain oap 
with radius around the point under considei'ation. In our particular case 
this can be done by subtracting from the measurements IT the mean value 
of the central block, that is for the inner zone, 


(1-19) 


N " Np - 0. 


The error function of the geoidal undulations in the inner zone is then of course 
likewise zero. On the other hand the neglected potential function itself will 
generate an error. 

We win give now a summary of the results of the error considerations. 
All the derivations in the further sections are done only for point values. Under 
some additional assumptions we will enlarge these results also to the error of 
mean values. The point error can be separated into three terms* 

(1-20) m^{Ag) = m^(Ag„^ ! m®(Ag„) -i- m^(Ag„) 


m (Ag) ........ mean point error 

m(Agn). mean error due to data error and the 


neglection of the inner zone. 


S'- 



m (Agsa) mean error due to unlcnown or 

variable sea surface topography 

m(Ago;t) mean error due to the neglection 

of the outer zone. 


Under these three types of individual errors, ra{Agsj) seems by far the most 
critical. Therefore, we will consider tliis error first. 

The following formulae and numerical values are dependent on the main 
assumption that mean geoidal undulations in blocks of size b[km] + b[km] and 
with an accuracy of m(N) [cm] arc given. The number 0.34 in the second 
term of (1-17) depends on an estimate of the root of the variance of the gravity 
gradient of 


11-^11 = i24EU, 


EU Eotvos unit 


drawn from (Metzger and Jircitano, 1974) 


Under these assumptions the accuracy of point anomalies at the midpoint of the 
blocks can be estimated by 


(1-21) m®(Ag„) = 


0. 34^ b^ 


(m(Agn) in [mgal], m(N) in [cm], and b in [km] ). 


The first term of the right hand side is due to measuring errors m(N), the 
second term is a result of' the neglected potential function In the inner zone. 
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With the aid of this formula we will now give aa estimate of the accuracy 
of mean anomalies of block size 

c > b/- 


In order to do this we make two additional assumptions ; 

1) The error of representation for sufficiently small blocks 
of size b is small as compared with the point error of 
the midpoint due to measuring errors, hi this way we 
can consider (1-21) also as an error estimate for the 
mean anomaly values of small blocks. 

2) The error values between neighboring blocks are only 
sightly correlated. If they are not correlated we get 
as error of a mean anomaly of block length c 

lh)(Aga) =-^m(Agn), 
o 

if they are very strong correlated we get 
mo(Ag„) = m(AgB). 

As an estimate of the error of mean anomaly values of block size c we obtain 
therefore 

(1-22) ^m(AgB) 5 nTa(Agn) < m(Ag„) . 

c 

The correlation between errors of neighboring blocks depends on the 
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correlation of the gravity gradient in the midpoints and on the correlation due 
to the measuring error. It is known that the covariance function of the gravity 
gradients goes very fast to zero (Moritz, 1974). According to the results of 
Gopalapillai (Gopalapillai, 1974, p.77), the error covariances between l®xl° 
blocks are also very small. So it seems not unlikely that the upper error limit 
in (1-22) can be reduced further. 

The following two examples give an idea about the accuracy of the 
mean anomalies. In order to get a satisfactory ratio.between the two error 
sources (that is, due to measuring errors and due to the neglected innor zone) 
we assume for both terms on the right hand side of (1-17) the same magnitude. 
This is attainable by a block length of 


(1-23) 

With 


b - 


I 


6, 1 • m (K) 
0.34 


m(N) = ±100 cm, b == 42 km 
we get according to (1-22) 

e = 1° X 1^ ; t8mgul < nT<,(Aga) < ±20mgal 

c - 2” X 2° ; ± -1m gal < mq(AgQ) < t:20mgal. 


In the case of 


m(N) = ±10 cm, b = 13 km 


we get 
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c 


; ilmgal ^ (^gn) ^ ±6mgal 


These results are in good agreement with the results of the error considerations 
in (Rapp, 1974). 

We will conclude this introduction with a few comments on the two other 
error sources. In section three it is assumed that the sea surface topography is 
known with the oKception of meandering currents such as the gulf stream. An es- 
timate of the error due to meandering currents is very difficult. A very rough 
estimate gave 


iSmgal £ m(AgBa) < i? mgal. 


These figures are point errors. However a systematic behavior, which would 
result in errors of the mean anomalies of the same magnitude, cannot be excluded. 
Hence attention is required in areas of strong boundary currents such as the gulf 
stream. 

If the gravity field of lower order is known the range of integration can 
be drastically reduced without remarkable additional error, ha this case it is 
possible to compute likewise gr-avity anomalies near to the coast. However, the 
remaining small errors are very strongly correlated at neighboring blocks. 

The error for several different circles \j)Q are given for a- gravity field 
known up to no = 2, no = 8 and no= 16 in section 6. In any case the influence of 
the constant term in the neglected part 

TT 2n 

Ag = - -~2s — r f M(0)(N ~ Np)sin0 d0 dry 

4rra : J,i J n 
tp-i/jo a=0 


may become fairly large and must be taken into accoimt in some way. 

The results given above are of com’se only a first estimate. They can- 
not replace test computations with true altimeter data and perhaps a comparison 
with accurate gravity anomalies, but they may give a good insight into the ditfi- 
cultities and the limits of the solution. 
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2« On the Solution of the Siverse Stokes' Problem in Pr<?>cise Gravimetric Geodesy 

The content of this section is entirely based on the results in (Lelgemann, 
1975), In that study the gravity anomalies at the earth's surface were reduced to 
gravity anomalies at the sphere with radius a. Then, the disturbing potential at 
this sphere was computed with Stokes' formula, combining the result with satellite 
derived information about the gravity field. From the disturbing potential at this 
sphere the disturbing potential on the earth's surface was computed with the quasi- 
geoidal imdulations, which can be found by the use of Bieiins equation 


S - — 
y 


C * quasi-geoidal undulation 

y normal gravity 


(In the present section we have to distinguish between the quasi-geoidal undulation 
t and the quantity defined by (2-19), The gooidal undulation N, used in all 
other sections, can then be Identified with N. or can be considered as C ia 
spherical approximation). 

In this report we are going precisely the reverse way. We transform 
quasi-geoidal undulations (as a result of altimeter measurements) into the dis- 
turbing potential, reduce this disturbing potential to the disturbing potential at the 
sphere a, solve the inverse problem of Stokes for the sphere and reduce the com- 
puted gravity anomalies to gravity anomalies at the surface of the earth. 


First we shall compile some well-known integral formulae. All these 
formulae are valid if the range of integration is a sphere. We consider a spherical 


polar coordinate system with r 
tude) and a fixed sphere with r 


(radius vector), 
- a = constant. 


0 (polar distance) and X (longi- 
On this sphere a continuous func- 


tion X should be given. An evaluation of x in a series of spherical harmonics 


should not contain terms of zero and first order, 
fimetion y on the sphere by 


Now we can define a second 


( 2 - 2 ) 


y 


a 

4it 


J X S(0) do. 
o 
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This is Stokes' well-known formula applied to a fictitious function x^ The 
function is given by 


(2-3) S(ib) = ir-r- 6sinJij/2 + 1 - Scosil) - 3cos^^>^7^(sin0/2 + sin®j^)/2) 

^ ^ ' sin(0/2) 

or as a series of spherical harmonics 
(2-4) S(t!J) = 2 

I) — 

We want a formula for the inverse problem, i. e. we will compute for a given hinc- 
tion y a function x. For this wo can consider the integral formula (2-1) as a 
Fredholm integral equation of the first kind. In a general case, the solution may 
become a very complicated problem, but in the case of the sphere wo get the solu- 
tion easily by a development of the functions involved in*o a series of eigenfunctions, 
that is, of spherical harmonics. Duo to the orthogonality relations of spherical 
harmonics we have 

CO m 

( 2 - 5 ) 2 y. = a ^ X. 


or by a comparison of the same terms of both sides of the equation 


(2-6) ax„ (n-l)y„. 


For the similar series 


(2-7) 2 

n = 3 




n 


yn 
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or 


(2-8) 


a * Yn 


a closed solution Is Icnown by the integral (Heiskanen-Moritz, 1967, p. 39) 


a'' 


2 n 


r da 


a 


(2-9) z - 

with 

i - 2asini/)/2, 

provided y is differentiable at least twice at the point P under consideration. 
We will write (2-9) in the following form 


( 2 - 10 ) 


with 


( 2 - 11 ) 


J M(ii!) (y - Yp) da 
a 


M(^^)) 


(2 sin {/)/2) 


TT • 


Writing equation (2-6) in the following form 


(ax + y)n - ny„ 
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we obtain 


or 


ax y 


4tt 


J M(i^) (y " Yp) do- 
er 


(2-12) 


X 


L- 

a 4Tia 


J M(ijl)) (y - Yp)da 
a 


or in series foi’m 

l33 

(2-13) X = ^ ^ (n-l)y. . 

Jl : ' n 


Assuming in this section that the results of altimetry are quasi-geoidal undulations 
C we neglect for the moment the influence of sea surface topography. We shall 
return to the more realistic case in the next section considering then the sea sur- 
face topography as an error source. 

As mentioned above the computation of Ag values can be managed in 
precisely the reverse way than the computation of geoidal undulations in (Lelgeniann, 
1975, section 4). We will use here the same notation and abbreviations. We solve 
the problem lilcewise in eight successive steps; 


1) Transformation from geoidal undulations C into 

the disturbing potential T at the surface of the ocean by 


(2-14) T - c • y 


y ..... . normal gravity 

formulae for y are given in (lAG, 1970, p. 58) 
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2) The correetion due to the indirect effect of atmospheric 
gravity reduction was computed to be maximal 0.6 cm. 
So we can set 


(2-15) 



3) Correction due to the. upward continuation term because 
of topography 


(2-16) 


Te - T - 
at^ - -H* Ag 


H ...... orthometric height 

Ag gravity anomaly 


This correction is necessary in the neighborhood of coasts and on islands. 
On the sea we have, of course, due to H-0 


6to = 0. 


4) Correction because of the ellipticity of the reference 
surface. The computation of the disturbing potential 
at the sphere with radius a is found from 

(2-17) T^ = T, - at3 


where 
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(2-18) 


COS® 9 • T (0, X) 


6t 


3 


e 


t2 


4 


5) We write in the following form 


(2-19) 


• 7« . 


7b normal gravity at equator 


In order to avoid the definition of a geoidal undulation in space, we will con- 
sider Nr as an auxiliary quantity defined by (2-19). However, Na can be 
considered as the quasi-geoidal undulation in spherical approximation and 
formula (2-20) may be more familiar then a similar formula, using <-ho poten- 
tial itself instead of Nr , 

(2-20) Ag, - . ( Nr + ^ J M(!/>) (Ni - (Nr)p) d(T). 

O' 

This is formula (2-12), applied to our special problem. 

6) Computation of the gravity anomalies at the ellipsoid 
from the gravity anomalies on the sphere with radius a; 


(2-21) Age - Agi -I- Bga . 

From (Lelgemann, 1975) we find 

a “ II 

(2-22) 6g^ ^ ^ ^ (n - 1) [C„t.R„„(0, X) H- DnaSnB(0, X)] 

n = d n = 0 

With 
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uimormalizGd spherical harmonics. 
They have the same definition as in 
(Heiskanen- Moritz, 1967) 


The coefficients C„q and Dnn can be computed in the following way. When 
(2-23) T(0,X) = OB Him (8,X) + H„,S„,(0,X)], 

n — H D=0 

then 

(2-24) Cna ■ A(n_gja Pjjb I' ■^(n+a)ti^no 

i^nn ~ H(n_a)a Pno ^ HnaQim H(o40)BiI*iiiti 

where 

(5n - 17) (n - m - 1) (n - m) 

4{n- l)(2n -3)(2n-l) 

- 6n^ - 8n^ + 25n + Gnm’^ + Gm^ + 21 

4(n ~ l)(2n i 3)(2n - 1) 

(5n -I- 11) fn ^ m I 2) (n -i* m + 1) 

4(n - 1) (2n + 5) (2n + 3) 


(2-25) PbB = 


Qtun 


rne 


7) The gravity correction due to topography is zero, because 
we want gravity on the sea surface. We hjive 
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(2-26) 


Ags = AgE + 


6gs 


with 


%a = Hp * Li(Ag), 


where Li (Ag) is an abbreviation, for 


Li(Ag) 


2tt 


f 


Ag - Aa, 


dcr. 


This correction is only different from zero, if we want to compute gravity 
anomalies at land. At sea, Hp Is zero and therefore we get always 


6ga = 0 . 


8) Add to this value the gravity effect 6g\ of the 
atmosphere 


(2-27) 


Ag - Ags + 6g^. 


Because this correction is constant at mean sea level we have 


6gi = -0. 87mgal. 

However the recommendation at the lAG (lAG, 1970) may be notified, that 
this correction may be handled as a measurement correction. In this case we 


20 - 



have of course 


6gi == 0. 

There are no difficulties in the computation of these correction terms. So wo 
can restrict pur accuracy considerations to the case of a spherical problem. 
Only the problem of the unknown sea surface topography needs a special con- 
sideration. 
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3. An Estimate of ^he Error Influence of Sea Surface Topography 


As Moritz (1974a) has shown we can compute the sea surface topography 
with geodetic measurements by a combination of altimeter data and gravity obser- 
vations. However, if we want to compute gravity from the altimeter observations 
the difference between the real ocean surface and an squipotential surface of the 
earth’s gravity potential must be known. 

The difference between the instantaneous sea surface and an equipotentlal 
surface of the gravity field is generated by many different forces. Therefore, the 
computation of this quantity may be very complicated (Wemelsfelder, 1970). For 
our present purposes it seems sufficient to consider the instantaneous sea surface 
as an equipotential surface disturbed by irregularities of four types ; 


1) very short periodic irregularities (e. g. ocean 
waves, swell) 

2) periodic or quasi-perlodic irregularities (e. g. tides) 

3) qua si stationery irregularities, which rclain their 
form but change their places (e. g. duo to the 
meandering Gulf stream) 

4) quasistationery Irregularities, which retain form 
and place (e,g. due to equatorial currents). 


From a geodetic point of viev? we consider sea surface topography only 
as an error source. So we will make two suppositions in view of the following 
error considerations; 


a) The very short periodic irregularities and the periodic 
and quasl-periodic effects by tides and seasonal effects 
are corrected either by statistical methods or by suitable 
models of the disturbing forces. 



b) Irregularities of type three and four must be given 
by oceanographic science or we have to consider 
Ihese irregularities as a possible error source. 


Because of the second supposition, the error discussion was made under 
the thesis: If oceanography does not need the help of altimetry for a solution of 
their problems the altimeter data can be corrected to geoidal undulations. If 
oceangTaphic sciences want altimeter data for a computation of certain phenomena 
we will consider these as an error source in our model. 

From this point of view the difficulties seem not so important with irrogn- 
laritles of t 3 ;pe four, which result from the shape and the magnitude of the larg'S 
scale- features of the mean sea surface. Information is already given in form of 
maps, from which we will mention the map from Stommel (Sturges, 1972) and the map 
from Xiisitzin (Bapp, 1974), The difference of about 1,5 m between the values in 
the two maps Is largely constant. An estimate of ±30 cm of the uneei^tainty of these 
maps is given by Stui*ges, He mentions that this error could probably be i*educed 
to approximately ± 10 cm. However, the difference between oceanographic and 
geodetic results in the large scale features of mean sea level along the coast lines 
are not yet clarified. 

The situation seems not so well defined in the case of strong boundary 
currents such as the gulf stream and maybe in the case of eddies. The tracking 
and also the detecting of strong boundary currents and of eddies could be one of 
the tasks which oceanographic science eicpeotB from altimetry (Sturges, 1972). 


Strong boundary currents such as the Gulf stream have the following form 
(Stommel, 1965), 
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The current flows perpendicular to the plane of the page. The height difference 
between the warm water of the sargasso sea and the colder water at the north 
side of the stream is ejqjected to be up to two meters. So we have to expect a 
slope of the stream of 


2 m: 100km = 2 • 10"® . 


The main problem is the meandering of the gulf stream. Although the meanders 
can remain in a fixed position for a short time, they usually change position sub- 
stantially during a time scale of perhaps two weeks. An expression of this mean* 
dering is seen in figure three in (Sturges, 1972), 

It seems very difficult to obtain an exact knowledge of this mean- 
dering, which is expected in the eastern part of the gulf stream to an horizontal 
extent of 500 kilometers. So we have to consider the sea surface topography 
generated from such a current as an data error fg . 

From formula (1-2) we get 

(Ag 4- 4) - -{^(N + fg) + J M(J/J) [(N + fg) 

a 

- (IS^ 4(^)p )] da 


or 


^ M(i|)) {%i - (Sg)p) da. 


(3-1) 


% 


a 


^ “ 


7a 

4na 


An error estimate based on this formula seems not easy. We will give a rough 
estimate of this error using in a completely different way our knowledge about 
the slope of the sea surface topography. 

For the ratio between the norm of the slope Av of the geoid and the norm 
of the gravity anomaUes Ag we draw from (Melssl, 1971* page 42), 
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y<, 

/r- 


( 3 “’ 2 ) 


£ 





where we have let G - y, . 


We will apply this formula to the error fj. We do not know the variance 
of fv , but only a maximal value of max(fr ) = 2 • 10 . We interpret the error 
due to the sea surface topography of the boundary currents such as the gulf stream 
as an stochastic process with a normal distribution on the sphere. Ih this case 
we have 


(3-3) max | fg 1 === 3 • |j jj = 3 • m (Aga») 


and 


(3-4) max I fv I = 3 • 1 fv 1| = 3 • m(Av«s). 


Then we get the following inequality 


..T"., 5nax I fv I < max j fg j max | fv 1 . 


With 


and 


max i fv I == 2 • 10”^ 

% ^ 980000 mgal 


we get as an rough estimate 


- 25 - 


8 mgal £ max [ f g { -e 20 mgal 


or 


(3-5) 


i 3 mgal 5 m (Agsg) ^ ±7 mgal . 


Because of the assumptions above the derivation of this estimate is 
very wealc, but it gives at least an expression of the magnitude of this error. A 
recovery of gravity data in areas of strong boundary currents should be done with 
uttermost caution, including a very carefully error analysis of data and results. 
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4. The Error Influence of the Neglected Inner Zone 


We have assumed already that our data are mean values of geoidal 
undulations N of blocks of size Ft, = b blkm®] with an root mean square 
error of m(N), We will compute the value Agp of the gravity anomaly always 
at the mid-point P of the blocks. For the value (N - 1^) inside of the inner 
most block we get 


(4-1) 


(N - Np) - (Np - Np) « 0, 


replacing in the integration procedure the function N by its mean value N. 

We will now estimate the error m (Agi) generated by this procedure, 
hi view of the fact that we make only an error consideration, the integration 
over the inner square block may be approximated by an integration over an cir- 
cular cap of the same area with sufficient accuracy. Because of 


(4-2) 


we get 


(4—3) 


Fi,= b?icb = Si^TT 



as a value of the radius Si of this circular cap. Approximating the radius H-t 
of a spherical cap by the radius s ^ of a plane circle 

(4-4) j&i - 2asin)|)i/2 = asin^i ^ Sj 

and the spherical surface element by 
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(4-5) 


d(T = sinijj dip do/ = s ds da, 


we get for the influence of the inner zone 
(4-6) Agi = - 


2 tt ipi 

— 1 f M(^) (N - Sin*/) di/j da 

^J,j 


2tt Si 

■i" I I (N - N.) ^ ds dor . 

a=0 s = 0 


We develop the geoldal undulation N at the point P into a Taylor series 


(4-7) 


1 / 


N - Np -h kN,; + yNy (x-^Nxx f- 2xyKxi- + y'^Nyy) U.. 


with 


Nx -(-^“)p 


and 




Introducing polar coordinates, 


X = s • cos a 


and 


y = s * sin a 


we get the following expression 


N = IS^ + s(NxCosa -h NySinO!) + (*-^)(NxxOOS^a 


+ 2 ISIxy cosa sina + Nyy sin®o) . 
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The integration over the azimuth Of gives 

Si 

Agl ^ f (N«X + NyyjdS 

^ S=0 

and the integration over the distance s 


Agj - - 4 ~ ^ '~ 


or 

(4-8) Agi = - 

Kt t? • • • ' 


% • -gi ( 

4 \ r^X ?>y 

. . deflectionB of the vertical 


We cannot estimate the error Influence from this formula, because 
we have no statistical information about the size of the horizontal gradient of 
the deflections of the vertical. However, based on Laplace's equation we can 
exchange this expression into 


(4-9) 


^ = -G (4^ + 

rir \ X 


?sy } 


(Heiskanen-Moritz, 1967, (2-221) ). 


Beeause 



we obtain for the influence of the inner zone; 


(4-10) 


Agi 



This expression is better suited for an estimate of the error which appears with 
the neglect of the inner zone* The variance of the gravity gradient is known to 
be (Metzger and Jircitano, 1974); 


II II = ±24 EU. 


EU ...... Eotvos unit 


With the relation 


1 Eotvos unit = 1 EU = = 10"® sec ® 

10 km 

we get as an estimate of the root mean square error of the point anomaly for the 
midpoint of the inner zone due to the neglecting of this zone 

24 1 

m(Agi) = ± • b [mgal] 


or 


nft(Agi) = ±0. 34 • b[mgal], 

where b must be inserted in the dimension km. It is easy to see that the value 


- 30 - 


0. 34 depends very strongly on the estimate of the variance of the gravi*'y 
gradient of ±24 EU. 

It may be pointed out that the formula (4-10) is valued only for a suf- 
ficiently small zone, dependent on the desired accuracy. Eor a more detailed 
examination see (Meissl, 1971a). 
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5a Error Cotisideratioas for a Given Accuracy of the Altimeter Data 


As in the last section we assume that our data are mean values of geoidal 
undulations N of blocks of the size Ft = b 4: b with a root mean square 

error of m(N), Furthermore we make the assumption that these measurements 
are homogeneous and uniformly distributed over the whole sphere. The more 
realistic case, that the mean values of the geoidal undiilations are given only in 
a certain neighborhood of the point under consideration, is discussed in the next 
section. 


We have to study the error propagation in the inverse Stokes ' formula 


(5-1) 


Ag 


^ J (N-Np)M((S)dcr. 

^ 0 


For the moment we shall consider only the integral part of formula (1-2), taking 
the linear term into account afterwards. 


The main difference as compared with Stokes’ formula is the stronger 
singularity of the kernel M0) at the point of interest. But we have already 
excluded a neighborhood around this point and the error estimated in the last 
section. So we can apply for the error considerations the same methods as for 
the error considerations in the case of the computation of geoidal undulations 
from gravity anomalies (Heiskanen-Moritz, 1967), The individual error ^ 
of Ag is 


(5-2) 


e = 


ya 

4Tia 


2 n 

O'io 



(j/), 0 ) M ($) sin0 di/j dor 


where p(0, a) is the individual error of the geoidal undulations. 

We consider only the error fimction l', talking the constant term 
into account with 
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7^ * 


2rr TT 

r r M (J/J) sin^^ di/) da = 

4,TTa qjJ,q ^1 q 


(See Heiskanen-Moritz, 1967s P* 37). Errors of fhis size can be neglected as 
discussed at the end of the section. 


The st^are of (5-2) becomes 

„ 2n TT 


“ f- ;, f f (i|j, a) M (</)) sin(/) diji da 

' O^o Olio 

2tT TT 

• I j a’) sln0’ dJ/j’ da’. 

fv4o 0'=O 


Because products of definite integrals may be written as one multiple integral 
we get 




-Is- 

4rra 


2n 

a=o 0 


TT 


) J I 


2tt 

I 


r P(i|), a) v(i/)', a’) 


• M(0) M(i^)') sinii) sin0’ d0 dad^’ da*. 


We form now the average of both sides of this equation arriving at 


2tT TT 


(5-3) 


m‘ 


2rr TT 
P 


'(^g) = (t^) f f f 1 Of, 0’, a’) 

^ a^ 0=0 a^=0 0*^0 


• M (0) M (0*) sin0 sin0' d0 da d0' da'. 
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In this formula, m^(Ag) is the error variance of the computed Ag-values and 


(5-4) ( 7 ( 1 ^), a, il3\ a') = or(6, X, 9% X’) = { v{6, X) v(0', X') } 


(the braces mean the average over the whole world) 


is the error covariance function of the geoidal imdulations. This is the general 
form of the error propagation in Stokes' inverse formula. Of course this formula 
is too complicated for our purpose to get an estimate of the order of the error. 

So we include two assumptions which can be drawn from the presuppositions about 
the measuring data; 

1) Only errors at neighboring poihts are correlated. 

We can make this assumption if errors of different 
blocks are uncorrelated. The error covariance 
function of the geoidal undulations may then assumed 
to be zero for points fiirther than b • /T apart. 

2) The accuracy is the same for every point of the 
earth surface. If we work exclusively with altime- 
ter observations this assumption is only valid at the 
oceans. For this reason the outer zones must be 
excluded from the computations. The additional 
error due to this exclusion will be considered in 
the next section. 


By as sump lion one the integrant is noticeably different from zero only 
if (0* = 6) and (X* = X), because the error covariance function may be considered 
as zero for two points more than h/T apart. Hence, we may approximate 


(5-5) 


g tt 2tt n 

m®(Ag) = J J J J O0ta,il)\OL') 

^ aio i/j=o 0'=O 


• M®(0) sin|/J siinl)' d</) do'di^i' da' 


and perform then the integration over J/ ' and n>' first. Wo define 
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(5-6) a® J J iji\ o:') sin0' di/j' dn:' = E(0,ry) 

«'=0 0*=O 

E , o . . . error constant of the altimeter dhta. . 

As a maximal limit of the real error covariance function, we will assume 
that the error covariance function is const = m(N)'' inside of each block and zero 
outside. 



as a maximal estimate of the error constant. . Because of assumption two the error 
constant E is independent of position. In this way we obtain 


(5-8) 


a 

m 



TT 

r [M0) sinij} dip da, 
ipio 





way. Consequently, we exclude the origin by beginning the integration, with 


0 = 01 

m^(Ag) = (- 

or 


(5-9) 

m^(Ag) = — 

with 

. n 

(5-10) 

Ji(0i) = J 
0 = 0 

Because of 



E 


8tt 




[M(l/i)]^ sint/i dS}i 


E 


8tt 


(-^J • Ji#i) 


MW = 


2^ sin® ij}/2 


we get 


n 


Ji(W = -lT,I,'sin>/2 


in order to solve the integral we introduce the new variable 

z - sinij[t/2 
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with 


sin tpi/2 = t 


and 


sini^i di/) = 4z dz. 


We ean easily verify 


(5-11) J.(M = ^ blATF- 


In this way we get as a final result 


(5-12) ni®(Ag) - — r t t - t : “ 1 1 

^ ^ ^ ^ 128tt \ a® / L sin^<^i/2 J 


We will now specialize this formula for small angles 0i, We get 


(5-13) 


m^(Ag) = 


E 


128 ri 




Sin* ipt/2 


Taking into account (4-3) we can approximate 


sin0i/2 = 0i/2 = St/2a = b/2/rT 


or 


16 TT® a* 


sin* tpf/2 
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Furthermore, we have jErom (5-7) 

E = m®(N) • . 

Putting these ejqDressions into (7-13) we arrive with 

= -5- • “"W 

or 

(5-14) m(Ag) - J m(N). 

This can be wri'-ten in a form 

(5-15) m(Ag) = • m(].‘q) 

( b in [kml, m(N) in [cm], m(Ag) in [mga]] ), 

Taking finally also ^he error due to the neglecting of the inner zone 
into account we get 

(5-16) tti®(^gB) = m®(Ag) + m^(Agi) 

or 

(5-17) m^(Ago) = - - m -f- 0.34® b'^ 
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(b in [km], m(lSl) in [cm], m(Ag 5 ) in [mgal] ), 


This is the variance of the error of the gravity anomaly at the midpoint of a 
compartment. 

We have neglected up to now the influence of the constant term in the 
inverse Stokes’ formula (1-2) 


(5-18) N = 0.15N[mgal]. 

a 

N has the dimension meter. It is easily seen that an error of 1 m 
in the geoidal undulation generates an error of 0.15 mgal in the result. This 
error can certainly be neglected in the error considerations. 
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6, ‘The Error Influence of Neglected Outer Zones 


In the last section it was pointed out that the derived error formula is 
only valid in the case of geoidal undulations given homogeneously over the whole 
earth. From altimetry alone there is of course a lack of data over the continents. 
In the present section we will estimate the error which occurs if we extend the in- 
tegral only over a certain area, a cup of radius 0 q, m this case we neglect the 
influence of the function itself outside of this cup. 

First we split up the integral into 


(6-1) 

Ag = Agi 

+ Ag(Np) + 

Ages 

with 








2tt 


(6-2) 

Agi = - 

Va r 

4 TT a J 

J 

M (^) (N - Np) sin(|l) d^ da 



lf)=0 01 

= 0 





TT 2tT 

(6-3) 

Ag(X^) = 

Np • — ^ 

4Tia 

0 = 

J M(^) sini^ dt^ da 
= 00 o;=0 



rr 

2tt 

(6-4) 

Agoi = 

- 

^ TT Cl .1 


[ M ((/))• N ’ sin0 d0 da . 



*± 1 1 cL , V 


a=0 


If we integrate only over the cup ip ij)o we neglect of course the 
terms Ag(Ii^) and Ag^j. First we will investigate the magnitude of Ag(Np). 
Starting the integration in formula (6-3) with the integration over a and 
inserting M(0) by the expression (1-3) we get 




Np 


n 

Vo 1 


sin^ dJ^) , 
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In order to solve this integral we insert 


sin il)/2 = z 


which lead to toe result 


(6-5) 


Ag(Np) « n. 


[ 


2a L sin i{)o/2 



The following table gives an Impression about the magnitude of this correction 
term for various cap sizes j/Jq. 


Table One . (Np in meters, Ag in milligal) 




2° 

4,3 

Np 

5° 

1.7 

Np 

10° 

0.8 

Np 

20° 

0.4 

Np 


[meter] geoidal undulation at point p 

Ag(li^) [mgal] error due to toe neglecting of toe 

constant Np 


Because this correction term may become fairly large it must always be taken 
into account in toe course of excluding the outer zones from toe integration. After 
this we will compute toe magnitude of toe neglected term Agoj given by 




y«. 

4Tra 


n 2 tt 
(y-0 


N • 


siu0 d<& dcy , 
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We in<:roduce disconHnuous funcf'ion 


(6-6) 


M(0) 


f 0 if 0 ^ 0 < !/Jq 

I M(iiJ) if 00 ^ ^ ^ 


getting 

(6-7) 


IT 


2tt 


^go 


4Tta 


J 


0 = 0 0=0 


N • M(0) sin0 d0 do, 


The function M(0) can be developed into a series of Legendre polynomials 


(6-8) M(0) = y — y ■- U, P„ (cos 0). 

^ A 

H“0 


For the moment we will assume the coefficients U„ are known. A formula for the 
computation of the U„ will be derived aftei’Wards. Inserting the series (6-8) into 
the integral (6-7) and interchange the order of summation and integration we get 


(6-9) 


Agdz 


% 

Stra 



n 


(2n -I- 1)U„ 


0 = 0 


2n 



N 


Pn (COB0) sin0 d0 dor. 


We can write N as a series of spherical harmonics 


CO 

N = y N,(0,x). 

n 


As as solution of the double integral we obtain 


TT 2n 

[ r N '• P„(cos0) sin0 d0 dor 

0 = 0 or=0 


(2n *i‘ 1) 
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and with this intermediate result 


( 6 - 10 ) 


= - 


2a 


CD 

y UnN„(0, X). 

o 




Expression (6-10) may be considered as the error at a certain point p = P(6f X), 
caused by neglecting the geoidal undulations beyond a circle of radius ipo around 
this point p. However, what we want is not an individual error but its variance 

(6-11) ||ig« ir = = -^ { (£u.N.(9,X)) }. 

»=a 

Starting with this expression we get, similar as in (Heiskanen- Moritz. 1967, p,261), 
after some manipulations 


2 CO ra 

Uso. 11" = V I UnUn- {N„N„.}.. 

ti=:r’ n’riS 


It is well-known that all {Nn NnO are zero except for the case n' = n. Therefore 
we obtain finally 


(6-12) 




ii = a 


By Kn we have denoted the degree variance coefficients of the disturbing potential 
respectively of the geoidal undulations itself. 


It remains the evaluation of a suited expression for the coefficients Un. 
We start with the integral expression 




4 


M (^ii) Pn (COS0) sin!^ dj;!) 
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or 


TT 

r 


P„(cos^) sin0 dlf). 


(6-13) 





00 


1 

sin^ 0/2 


The main difference to the similar problem of the computation of Molodenskii's 
coefficients Q„ (Heiskanen-Moritz, 1967, p. 260) consists in the fact, that the 
kernel M(0) is much easier than Stokes' function. Introducing the new variable 

z = sin 0/2, t = sin0o/2 


COS0 = (1 - 2z®), sin0 d0 = 4zdz 

we obtain, 

t 

(6-14) U. = - I P„(l - z=) 


We will develop the Legendre polynomials into a power series of z. For this 
purpose we start with the well-known series representation of the Legendre 
polynomials 


(6-15) 


P„(l - 2z^) = 2"" y (-1)*' 
k = o 


(2n - 2k) I 

k I (n - k) ! (n - 2k) I 


(1 - 2z^)“-2»' , 


r ^ INT (n + 2). 


Because of the binomial expansion 
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(6-16) (a-b)“ = y (-1) 


V=o 


we get 


n-gk 


(1 - 2z®)’^2k = y — (n - 2k) ! 3 gy 

' ^ ^ ^ iv! (n - 2k - U)\ 


V -0 


or 


(6-17) 


(1 - 2z®)""3>' 

„2 


n-Kk 


L ' ’ v\(n-Zk-V)\ • 


V =D 


With this transformation of Legendre polynomials we can write the integral (6-14) 
in the following form 


U„ 


ti-3k 


1 ksiO 




where 


•^nk 


(2n-2k)l 2-” 

^ ^ k ! (n - k) ! (n - 2k) ! 


and 


BnkV 


(- 1 ) 


(n - 2k) ! 2^ 
t (n - 2k - i^) ! * 


We interchange the order of integration and summation and because of 2(i^ - 1) ?^-l 
we can perform the integration with the result 
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( 6 - 18 ) u„ = - y Am, y 

u U 2[^-l 


it = o V=o 


We will give this power series a more convenient form. First we can write 


n-2tt 


Un = - y y (-if'’"*' — ^ [2n " 2k) ! 1 rf.(sv-i)_ ^-i 

, ^ ^ V ! k ! (n-k) ! (n-2k-V) ! 2u-l ^ ^ 


'.c = o v = o 


Here we have of course r = INT(n/2). Now we will change the order of sum- 
mation, that is, we will sum first over u. Because of the case k - 0, the sum- 
mation over V has the range from 0 to n. Because of f - n - 2k we have for a 
summation over k a range from k - 0 to k = p with p = INT( (n-v)/2). We arrive 
therefore at the desired result with the expression 


(6-19) 



v=o 


^ 1 

k=0 


(2n - 2k) ! 

l>l k ! (n-k) ! (n-2k-P) ! 


2p-1 '■ 


with 


P===INT( (n-i^)/2) 
t = sin0o/2 . 


The final result (6-19) was cheeked by an analytical integration of formula (G-14) 
up to an order of n = 3, 

For an estimate of the error due to the neglecting of the outer zones we 
use the degree variances from the covariance function (1-8), computed by Tscherning 
and Rapp 


Kn 


a® \ 425,28 

) (n - 1) (n - 2) (n + 24) 


5 


n+X 


with 


S ^ 0.999617, 
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Ip this case we get as error 


( 6 - 20 ) 


111 (Ag„) = 


425.28 


00 

y 




.a+i 


(n - 1) (n - 2) ()i + 24) 


[mgair 


A computer program was written in order to get numerical results. These results 
are given in the following table for various radius 0o of the outer zone and different 
degree Dq of a reference field. The dimension of m(Ago:.) is milligal. 


Table Two . Influence of the Zone Beyond the Radius $0 for Reference Model 
of Degree no.- 


00 

2 

8 

l(i 

2 ° 

10.2 

0.5 

0.1 

5 ° 

3.3 

0.1 

(.0 

10 ° 

1.1 

0.0 

c.o 

20 ° 

0.2 

0.0 

O.v 


If we imow exactly a reference field up to a degree nn = 16, then seems 
generally sufficient to carry out the integration up to a circle of - 2° . So 
it may be also possible to compute gravity anomalies near the coast of the con- 
tinents from the altimeter data. 
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7, Smnmary 


A detailed description of a model for the recovery of gravity anomalies 
from altimeter data with 10 cm accuracy is given in Section Two, The disturbing 
potential at the sphere with radius a is computed from the geoidal undulations at 
the earth’s surface. After this, gravity anomalies at this sphere are recovered and 
then these anomalies reduced to values at the earth’s surface. The practical solu- 
tion consists of small correction terms to the inverse Stokes’ formula. 

The influence of the sea surface topography is discussed especiaEy in the 
case of meandering currents such as the Gulf Stream. From the influence of the 
unlcnown meandering of these currents a rough error estimate indicates an expected 
error of 


±3mgal 5 m(AgaB) ^ i7mgal, 

m(Agafl) error of the gravity anomaly due to 

sea surface tppog-raphy 


The most critical error Influence is probably due to the data errors and due to the 
neglect ion of the inner zone. It was assumed that mean geoidal undulations are 
given all over the earth in blocks of size b and with an accuracy of m(N) , 

Point errors at the midpoint of the blocks can be estimated with 




6.1^ m^(N) 
b® 


+ 0.014® 



b® 


m(AgJ error of the gravity anomaly in [mgall 

m(N) error of the mean geoidal undulations in [cm] 

b block size of the mean geoidal undulations in pem] 
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mean value of Ihe vertical gravity gradient 
in Eotvos units 


For m{N) = ±10 cm, b == 13km and ^ j| = 24 EU we get an estimate of 

r 

m(Ag^ = ±6nigal. 

An estimate of the accuracy of mean anomalies of block size c, 

c > b , 


gives under certain natural assumptions 


m(Agj} £ me (Ag^ ^ in(AgB)‘ 


voc (Ago) 


error of the mean anomaly of 
block size c. 


With the numbers above and a block size c == 1'' we obtain 


±lmgal ^ me (Ago) ^ ±6mgal. 


The results are in good agreement with the results of the error considerations in 
(Rapp, 1974), which are based on least squares collocation. 

In the last section the error of neglecting an outer zone in the inverse Stokes 
integral is investigated. Such a procedure may be necessai’y because it is not yet pos 
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sible to compute geoidal undulations from altimeter data on land. If the harmonics 
of lower order up to a degree Uq = 16 are known, the range of integration can be 
restricted to a few degrees without a significant additional error. However, a cer- 
tain correction has to be talcen into account in this case. 
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